Long-Term Evolution of Massive Black Hole Binaries. II. 
Binary Evolution in Low-Density Galaxies 
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ABSTRACT 

We use direct-summation TV-body integrations to follow the evolution of binary black holes 
at the centers of galaxy models with large, constant-density cores. Particle numbers as large 
as 0.4 x 10 6 are considered. The results are compared with the predictions of loss-cone theory, 
under the assumption that the supply of stars to the binary is limited by the rate at which they 
can be scattered into the binary's influence sphere by gravitational encounters. The agreement 
between theory and simulation is quite good; in particular, we are able to quantitatively explain 
the observed dependence of binary hardening rate on N. We do not verify the recent claim of 
Chatterjee, Hernquist & Loeb (2003) that the hardening rate of the binary stabilizes when N 
exceeds a particular value, or that Brownian wandering of the binary has a significant effect on 
its evolution. When scaled to real galaxies, our results suggest that massive black hole binaries 
in gas-poor nuclei would be unlikely to reach gravitational-wave coalescence in a Hubble time. 

Subject headings: black hole physics — galaxies: nuclei - stellar dynamics 



1. Introduction 

Binary supermassive black holes are inevitable 
by-products of galaxy mergers, and their coales- 
cence is potentially the strongest source of gravita- 
tional waves in the universe (Thorne & Braginskii 
1976). Following their initial formation at sepa- 
rations of roughly a parsec, massive binaries are 
expected to "harden" as the two black holes ex- 
change angular momentum with stars or gas in the 
host galaxy's nucleus. The binary separation must 
decrease by 1-2 orders of magnitude if the black 
holes are to come close enough together that grav- 
itational wave emission can induce coalescence in 
a Hubble time. Whether this happens in most 
galaxies, or whether uncoalesccd binaries are the 
norm, is currently an unanswered question (Mer- 
ritt & Milosavljevic 2004). 

This paper is the second in a series investigat- 
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ing the long-term evolution of binary black holes 
in galactic nuclei. As in Paper I (Milosavljevic & 
Merritt 2003), we restrict our attention to galaxies 
without gas. In such an environment, a massive 
binary shrinks as passing stars extract energy and 
angular momentum from it via the gravitational 
slingshot (Saslaw et al. 1974). Early treatments of 
binary evolution (Baranov 1984; Mikkola & Val- 
tonen 1992; Rajagopal & Romani 1995; Quinlan 
1996) represented the galaxy as fixed in its prop- 
erties as the binary evolved, and inferred the bi- 
nary's evolution from rate coefficients derived via 
three-body scattering experiments. This approx- 
imation is reasonable during the binary's initial 
evolution, but once the separation has decreased 
by a factor of order unity, the assumption of a 
fixed background is no longer valid. The binary 
quickly (in a galaxy crossing time) interacts with 
and ejects most of the stars on intersecting or- 
bits, and any subsequent binary-star interactions 
require a repopulation of the orbits in the binary's 
"loss cone." This argument has motivated various 
hybrid approaches to binary evolution, in which a 
model for loss-cone repopulation is coupled with 
rate coefficients derived from scattering experi- 
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mcnts in a fixed background (Zier & Biermann 
2001; Yu 2002; Milosavljevic & Merritt 2003; Mer- 
ritt & Poon 2004; Merritt & Wang 2005). 

Two regimes of loss-cone repopulation can be 
identified (Paper I). If the time scale for encoun- 
ters (or some other process) to drive stars into the 
binary is short compared with orbital periods, the 
binary's loss cone will remain nearly "full," and 
the rate of supply of stars will hardly be affected 
by their loss due to ejections. On the other hand, 
if time scales for loss-cone repopulation are long 
compared with orbital periods, the binary's loss 
cone will be nearly "empty," and the binary's evo- 
lution will be limited by the rate at which new 
stars can diffuse onto loss-cone orbits. If the dif- 
fusion is driven by star-star gravitational encoun- 
ters, the binary's hardening rate in this regime 
will scale approximately as iV _1 (in a galaxy with 
fixed mass and radius) . For values of N character- 
istic of real galaxies, N J> 10 10 , encounters would 
be rare enough that the loss cone of a massive 
binary would remain nearly empty, and the hard- 
ening rate would correspondingly be very low. In 
effect, the decay would stall. 

While the hybrid models are informative, a fully 
self-consistent, TV-body approach to the evolution 
of massive binaries is clearly desirable. Scatter- 
ing experiments in a homogeneous background 
do not faithfully reproduce the interactions that 
take place at the bottom of a galactic poten- 
tial well, where a given star may interact more 
than once with the binary (Kandrup et al. 2003; 
Milosavljevic & Merritt 2003). If the goal is to 
follow the evolution starting from the pre-merger 
phase, when the two black holes were widely sepa- 
rated, A^-body techniques are unavoidable. How- 
ever most iV-body simulations of binary evolution 
published to date (Ebisuzaki et al. 1991; Makino 
et al. 1993; Governato et al. 1994; Milosavljevic 
& Merritt 2001) have been based on such small 
particle numbers that the binary's loss cone was 
kept essentially full by star-star scattering or by 
random motion of the binary. As a consequence, 
these simulations failed to reproduce the diffusive 
loss cone repopulation that would characterize bi- 
nary evolution in the large-A" limit and they can 
not easily be scaled to real galaxies. 

In this paper, we consider the A^-body evolu- 
tion of a massive binary in a very idealized galaxy 
model, with a Plummer (1911) density profile. 



The Plummer model has a large, constant-density 
core, rather different from the power-law nuclei of 
most galaxies. There are several reasons for this 
unphysical choice. (1) The low central concentra- 
tion of the Plummer model implies a long star-star 
relaxation time, hence our A"-body models are able 
to maintain an empty loss cone with fewer parti- 
cles than would be required if we had used a more 
realistic galaxy model. This allows us to approach 
more closely than heretofore to the diffusive loss- 
cone repopulation regime. (2) Also because of 
its low central concentration, the Plummer model 
evolves only slightly due to the influence of the bi- 
nary. This makes it easier to compare our A^-body 
results with a model in which the gross properties 
of the galaxy are assumed to remain fixed with 
time. (3) Our initial conditions are precisely the 
same as those adopted by Chatterjee, Hernquist 
& Loeb (2003) in their numerical study of binary 
evolution. These authors used a hybrid A^-body 
code in order to achieve large particle numbers 
and reached a striking, counter-intuitive conclu- 
sion about the A^-dependence of the binary hard- 
ening rate. Based on this result, Chatterjee et 
al. concluded that "a substantial fraction of all 
massive binaries in galaxies can coalesce within 
a Hubble time." As discussed below, we fail to 
confirm their result with our higher-accuracy in- 
tegrations, and reach a different conclusion about 
the likelihood of coalescence. 

The A^-body integrations presented here were 
the first to be carried out on a new, special- 
purpose computer that couples GRAPE hardware 
with a parallel architecture. §2 describes the com- 
puter, the A^-body algorithm and the tests which 
we carried out to verify its accuracy. The evolution 
of the binary is described in §3, and in §4 we show 
how the A^-body evolution can be reconciled with 
the predictions of collisional loss-cone theory. Our 
results are compared with those of earlier studies 
in §5, and §6 briefly discusses some implications 
of our results for binary evolution in real galaxies. 
§7 sums up. 

2. Models and Methods 

Our initial galaxy model was a Plummer (1911) 
sphere, with mass density and gravitational poten- 
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tial given by 



p(r) 



An rl 



(la) 



*(r) - GM gal (l + x 2 ) 1/2 ,x = r/r .(lb) 

Here M ga i is the total galaxy mass, ro is the 
Plummer scale length and G is the gravitational 
constant. We henceforth adopt standard iV-body 
units (Heggie & Mathieu 1986), G = M gal = 1, 
E = 1/4 with E the total (binding) energy; in 
these units, r — 37r/16 ~ 0.589. Particle posi- 
tions and velocities were generated in the usual 
way from the equilibrium, energy-dependent dis- 
tribution function f(E). 

Our galaxy model was the same one adopted by 
Chatterjee et al. (2003) in their A-body study of 
binary evolution. We also followed their prescrip- 
tion for introducing the massive binary into the 
galaxy: two point masses of equal mass, Mi = 
M 2 = M/2, were placed on nearly circular or- 
bits at r and — r, with r = 0.3. No adjustments 
were made in the stellar positions and velocities 
when introducing the black hole particles. We 
chose two values for the masses of the black holes: 
Mi = M 2 = 0.02 and Mi = M 2 = 0.005. Inte- 
grations were continued until a time of t = 250, 
corresponding to roughly 88 crossing times, where 
T cr = (2\E\)- 3 / 2 « 2.83. We carried out inte- 
grations with a range of particle numbers, TV = 
(0.05,0.1,0.2,0.4) x 10 6 , in order to test the de- 
pendence of the results on N. 

The sphere of influence of a (single) black hole 
of mass M is r h = GM/ct(0) 2 with a(0) the 
ID stellar velocity dispersion at the center of 
the galaxy. In our Plummer spheres, <r(0) = 
(2/3)^27^ « 0.532, yielding r h = (0.141,0.0353) 
for M = Mi + M 2 = (0.04,0.01). The scmima- 
jor axis at which a binary first becomes "hard" 
is given approximately by at — GM 2 /4cr 2 with 
M 2 the mass of the small component (Quinlan 
1996); for an equal-mass binary, an — rh/8 = 
(1.76 x 10~ 2 ,4.42 x 10~ 3 ) for M = (0.04,0.01). 
Thus, the two black hole particles moved initially 
on widely separated, nearly independent orbits. 

The A-body integrations were carried out on 
gravitySimulator, 1 a special-purpose computer 
cluster recently installed at the Rochester Insti- 



tute of Technology. This cluster contains 32 dual- 
Xeon nodes running at 3.0 GHz. Each node hosts 
a single GRAPE-6A accelerator board (Fukushigc 
et al. 2005) which can store data for up to 131,072 
particles and calculate at a speed of 125 GFlops, 
allowing the entire cluster to carry out simulations 
with particle numbers of ~ 4 x 10 6 and at speeds 
of approximately 4 TFlops. Communication be- 
tween the GRAPE boards and their host nodes is 
via the standard PCI (32bit/33MHz) interface; the 
PC nodes are connected via a high-speed Infini- 
band network switch with bandwidth of 10 Gbit/s 
(duplex) . 

The TV-body code was an adaptation of Aarseth's 
NB0DY1 (Aarseth 1999) to the GRAPE cluster. 
The gravitational force acting on particle i is 
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m k (n - r fc ) 
l (e 2 + \v l -r k \>) 
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where m» and are the mass and position of 
the ith particle and e is the softening length; the 
gravitational force constant has been set to one. 
The integration of particle orbits was based on 
the fourth-order Hermite scheme as described by 
Makino and Aarseth (1992). We adopted their for- 
mula for computing the time-step of an individual 
particle i at time t, 



At, = 



\ a (t)\\ a (D(t)\ + \k(t)\ 2 
a(t)||a(3)(t)| + |a( 2 )(t)| 



(3) 



Here a is the acceleration of the ith particle, the 
superscript (j) denotes the j'th time derivative, 
and 77 is a dimensionless constant that sets the 
accuracy of the integrator. Since the particle po- 
sitions rfc must be up-to-date in equation (2) they 
are predicted using a low order polynomial. This 
prediction takes less time if groups of particles re- 
quest a new force calculation at large time inter- 
vals, rather than if single particles request it in 
small time intervals. For this reason, an integer n 
is chosen such that 



n-l 



(4) 



1 http://www.cs.rit.edu/~grapecluster/ 



with Ati given by equation (3). The individual 
time step is replaced by a block time step AbU, 
where 

&bU =(l) ■ (5) 
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The parallel algorithm begins by distributing 
the A particles randomly and evenly between the 
p GRAPEs. At each time step, the active particles 
on each node (i.e. the particles whose positions are 
due for an update) are identified and their coor- 
dinates broadcast to all of the other nodes, where 
the partial forces are computed. The partial forces 
are summed on the head node where the positions 
and velocities of the active particles are advanced, 
and their coordinates are then updated in the in- 
dividual GRAPE memories. More details of the 
parallel algorithm, including the results of exten- 
sive performance tests, are given in an upcoming 
paper (Harfst et al. 2005) . Most of the integrations 
used eight nodes; integrations with the largest par- 
ticle numbers (A = 0.4 x 10 6 ) used either 16 or 32 
nodes. 

The A-body code has two important parame- 
ters that affect the accuracy and efficiency of the 
integrations: the softening length e and the time 
step parameter rj. Figure 1 shows the effects of 
varying e and rj on the evolution of a binary with 
mass M\ = M2 = 0.02 in a Plummer-model galaxy 
with A — 0.05 x 10 6 . The initial positions and ve- 
locities of the two massive particles were the same 
as in the "production" runs. The overall accu- 
racy of the calculation as measured by changes in 
energy \AE/E\ (Fig. 1(a)) depends most strongly 
on the time-step parameter 77; for rj = 0.01, the 
fractional change in E is a few percent or less. 
Figure 1(b) shows that most of the error in the 
total energy comes from the binary, which is very 
poorly integrated when 77 is as large as 0.03. This 
figure suggests that a value 77 = 0.02 or smaller 
is required for accurate long-term integration of 
the binary. The evolution of the binary's eccen- 
tricity (Fig. 1(c)) is also poorly reproduced when 
f] = 0.03. Changing the softening length e appears 
to have much less effect on the evolution of a or 
e, although of course e must be small compared 
with the smallest separation attained by the bi- 
nary. Based on these tests, we adopted r\ = 0.01 
and e = 1 x 10~ 4 for the final integrations. 

3. Evolution of the Binary 

Figure 2 shows the evolution of the binary's 
semi-major axis in each of the eight integrations. 
The A-dependence of the hardening rate is clear. 
At late times, when a <C a/i, the binary semi- 
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Fig. 1. — Dependence of the global energy con- 
servation (a), binary semi-major axis (b), and bi- 
nary eccentricity (c) on the A-body parameters rj 
and e, for integrations with Mi = M2 = 0.02 and 
A = 0.05 x 10 6 . Black lines: e = 1 x 10~ 5 ; blue 
lines: e = 5 x 10~ 5 ; red lines: e = 1 x 10~ 4 . 
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major axis obeys a(t)~ x ~ C±(N) + C2(N)t, i.e. 
the binary's binding energy increases almost lin- 
early with time, and the hardening rate is a mono- 
tonically decreasing function of N. An approxi- 
mately linear dependence of 1/a on time is char- 
acteristic of both the "full loss cone" (small N) 
and "empty loss cone" (large N) regimes (Paper 
I). If the loss cone were completely full, however, 
the hardening rate would be independent of N, 
clearly a poor description of Figure 2. 

We define the instantaneous hardening rate to 

be 

««-!(:)■ < 6 » 

Figure 3 shows s(i) for each of the integrations, 
computed by fitting smoothing splines to a~ 1 (i) 
and differentiating. The constancy of s at late 
times (t J> 150) is apparent; aside from some wig- 
gles, the dependence of 1/a on time is well ap- 
proximated as a linear function. To a good ap- 
proximation, we can therefore identify a unique, 
late-time hardening rate s(N, M) with each of the 
integrations. We computed s by fitting a straight 
line to a" 1 (t) in the time intervals 150 < t < 250. 
Figure 4 shows the results, plotted versus parti- 
cle number N. The TV-dependence of the mean 
hardening rate is approximately 

log 10 s w 4.5 - 0.81 log 10 N, Mi = M 2 = 0.02, 
« 2.7 - 0.33 log 10 N, Mi = M 2 = 0.00, 

or 

— ( - ) « 3.3 x 10 4 /V" - 81 , Mi = Mi = 0.02, 
at \a J 

« 5.0 x 10 2 /V" - 33 , Mi = M 2 = 0.005. 

The very nearly linear increase of the binary's 
binding energy with time seen in all these simu- 
lations was something of a surprise. While a~ 1 (i) 
is predicted to be approximately linear in both the 
full- and empty-loss-cone regimes, at least at times 
before the binary has removed much mass from the 
core (Paper I) , the theory on which this prediction 
is based is only approximate. It will be interesting 
to see if the constant rate of hardening observed 
here is a very general feature of binary evolution. 

Also of interest is the evolution of the binary's 
eccentricity (Figure 5). The two massive particles 
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Fig. 2. — Evolution of binary semi-major axis. 
Solid lines: Ml = M 2 = 0.02; dotted lines: 
Mi = M 2 = 0.005. 
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Fig. 3. — Evolution of binary hardening rate. 
Color scheme is the same as in Figure 2. Solid 
lines: Mi = M 2 = 0.02; dotted lines: M x = M 2 = 
0.005. 
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Fig. 4. — Mean binary hardening rates in the in- 
terval 150 < t < 250 as a function of N. Open 
circles: Mi = M 2 = 0.02; filled circles: Mi = 
Mi = 0.005. The lines are least-squares fit to the 
iV-body hardening rates. Crosses show the pre- 
dictions from loss-cone theory, as discussed in the 
text. 
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Fig. 5. — Evolution of binary eccentricity. Color 
scheme is the same as in Figure 2. Solid lines: 
Mi = Mi = 0.02; dotted lines: Mi = M 2 = 0.005. 



are introduced into the galaxy model on approx- 
imately circular orbits with respect to the galaxy 
center, but encounters with stars induce a non- 
zero eccentricity even before the time that the sep- 
aration has fallen to ah, and the eccentricity con- 
tinues to evolve as the binary ejects stars. There 
is an iV-dependence here as well, in the sense that 
e tends to evolve less with increasing N, although 
the trend is obscured by noise. In addition, for 
larger TV, the "initial" eccentricity (i.e. the value 
of e when the binary first becomes hard) tends 
to be smaller, due presumably to the smaller size 
of random perturbations from passing stars. Fig- 
ure 5 suggests that the eccentricity evolution of a 
binary in a real galaxy with much larger N would 
be very small, although larger particle numbers 
will be required to verify this conclusion. 

The center of mass of the binary wanders quasi- 
randomly due to encounters with stars, both dis- 
tant, elastic encounters (Chatterjee et al. 2002; 
Laun & Merritt 2004) and "super-elastic" en- 
counters in which the binary's binding energy is 
transformed into linear momentum during ejec- 
tions (Merritt 2001). Figure 6 shows this gravi- 
tational Brownian motion in the four integrations 
with largest N. These plots show the motion of the 
binary with respect to a fixed (inertial) frame; be- 
cause the iV-body model as a whole drifts in space, 
there is also a systematic drift of the binary's mean 
position. We attempted to "take out" this drift 
by computing the position of the binary with re- 
spect to the galaxy's density center at each time 
step. However the structure of the galaxy models, 
with their large, constant-density cores, made this 
difficult since the position of the estimated den- 
sity center tended to vary from time step to time 
step with an amplitude almost as great as that of 
the Brownian motion. In any case, one can esti- 
mate the amplitude of the "random" component 
of the motion from Figure 6, and it is quite similar 
to what has been found in other iV-body simula- 
tions with similar initial conditions (Chatterjee et 
al. 2003; Makino & Funato 2004), as well as with 
the expectations from encounter theory (Merritt 
2004) . Below we address the question of whether 
the binary's Brownian motion might influence its 
hardening rate. 

The presence of the binary affects the central 
properties of the galaxy, although only slightly. 
Figure 7 shows the evolution of the central density 
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Fig. 6. — Motion of the binary center of mass in 
the four integrations with largest TV. The posi- 
tion of the binary is measured with respect to the 
fixed (inertial) frame; the systematic change that 
is apparent in some of the plots is due to an over- 
all drift of the TV-body model. Black lines: xcm i 
blue lines: ycM] re d lines: zcm- 
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Fig. 7. — Evolution of the central density of 
the galaxy, defined as the mean density within a 
sphere of radius 0.2 centered on the binary, in the 
two integrations with TV = 0.4 x 10 6 . Small dots: 
Mi = Mi = 0.02; open circles: Mi = M 2 = 0.005. 
Dashed line shows the mean density within r = 0.2 
in the initial model. 



in the two integrations with TV = 0.4 x 10 6 ; here 
the density was defined as the mass in a sphere 
of radius 0.2 centered on the binary divided by 
the volume of the sphere. The introduction of the 
binary into the model at t = causes an impul- 
sive change in the local density, particularly in the 
case of the more massive binary, but the effect is 
transient and thereafter the density changes only 
by ~ 15% over the course of the integration, due 
presumably to ejection of stars. The evolution in 
central density was almost the same for these two 
integrations and depended only weakly on TV as 
well. Such modest changes in the central proper- 
ties of the galaxy models could hardly be responsi- 
ble for the strong dependence of binary hardening 
rate on M and TV seen in Figures 2-4. 

4. A Model for Binary Evolution 

In an infinite, homogeneous background with 
fixed properties (density, velocity dispersion etc.), 
the binary's hardening rate is given by 

.-!(!) -*Se (9 ) 

at \a ) a 

with p and a the stellar density and ID velocity 
dispersion; H is a dimensionless decay rate that 
depends on the binary semi-major axis, as well as 
on other binary parameters such as mass ratio and 
eccentricity. Scattering experiments give if as 16 
for a hard (a -C ah), equal-mass binary (Hills 1983; 
Mikkola & Valtonen 1992; Quinlan 1996; Merritt 
2001). If we set p and a to their values in our 
Plummer-model galaxy at t = 0, equation 9 with 
H = 16 predicts s ~ 35, independent of M and 
TV. By comparison, the largest, late-time decay 
rate found in these integrations (for M = 0.01 and 
TV = 0.05 x 10 6 ) was s « 14 (Figures 3, 4), not 
terribly different. However, the strong dependence 
of the decay rate on M and TV in the simulations 
is clearly inconsistent with this simple model. 

As discussed in Paper I, there is a natural ex- 
planation for the TV-dependence of the hardening 
rate. After roughly a single galaxy crossing time, 
the binary has ejected all of the stars on intersect- 
ing orbits. Subsequent interactions between the 
binary and the stars occur at a rate that is lim- 
ited by how quickly new stars can be scattered 
into the binary's loss cone. The latter is defined 
as the set of orbits with pericenters r p < Ka(t), 
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with K a number of order unity; this expression 
reflects the expectation that stars will be ejected 
from the binary whenever they pass a distance 
~ a from its center of mass. But the scatter- 
ing of stars onto low-angular-momentum orbits de- 
pends on the two-body relaxation time, hence on 
the mean stellar mass, hence on N, and so one 
expects the hardening rate of the binary to be N- 
dependent as well. 

As a first step toward testing this model, we 
asked whether star-star encounters might be so 
frequent in the A-body models that the binary's 
loss cone is maintained in a continuously popu- 
lated state. The standard measure of the loss-cone 
refilling rate in a spherical, steady-state galaxy is 
q{E), defined as 



J lc 



(10) 



(Lightman & Shapiro 1977). Here 5 J is the change 
over one radial period in the angular momentum 
of a star on a low-J orbit, and Ji c is the angular 
momentum of an orbit at the edge of the loss cone. 
A value q{E) ;§> 1 implies that the loss cone orbits 
at energy E are re-populated at a much higher rate 
than they are de-populated by the central sink (in 
this case, the binary black hole), and the loss cone 
remains nearly full. A value g C 1 implies that 
the loss cone is essentially empty, and repopulation 
must take place diffusively, as stars scatter in from 

J Jlc- 

We computed q(E) in our Plummer-model 
galaxies, ignoring any effect of the binary on the 
galaxy. We used the precise, orbit-averaged defini- 
tion of q(E) given by equation (8) in Paper I. The 
results are shown in Figure 8. The binary param- 
eters appear implicitly in g, since Jf c = 2GMr t 
with r t = Ka the radius of the capture sphere. 
Figure 8 shows q as a function of E for the eight 
different combinations of (A, M) in our A-body 
integrations. We considered two values for J; c , 
Jf c = 2GMa h and Jf c = 2GMa h /10; the first 
value is the approximate "size" of the loss cone 
when the hard binary first forms, and the latter 
value is appropriate when a has decreased by a 
factor of 10, roughly the maximum amount of de- 
cay seen in these integrations (Fig. 2). Figure 8 
suggests that the integrations with the higher- 
mass binary, M = 0.04, took place mostly in the 
"empty loss cone" regime, q 1, while the inte- 
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Fig. 8. — The quantity q(E) that measures the 
degree to which the binary's loss cone is repopu- 
lated by encounters. Color scheme is the same as 
in Fig. 2. Solid lines: Mi = M 2 = 0.02; dotted 
lines: Mi = M 2 = 0.005. The thicker curves in 
each group are for r t — , the approximate radius 
of the loss sphere when the binary first forms, and 
the thinner curves are for r t = a^/10, correspond- 
ing to an evolved binary. These curves suggest 
that the evolution of the more massive binary in 
our A^-body integrations took place mostly in the 
"empty loss cone" regime, while the evolution of 
the less massive binary took place mostly in the 
"full loss cone" regime. 
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grations with M = 0.01 were mostly in the "full 
loss cone" regime, q ^ 1. 

This result provides a natural explanation for 
the different TV-dependence of the decay rates 
in the two sets of integrations. When the loss 
cone is empty, q <C 1, re-population takes place 
diffusively, as stars are scattered by other stars 
onto low-angular-momentum orbits. The scatter- 
ing rate scales inversely with the relaxtion time, 
or as <~ 1/N. Indeed, the binary decay rate for 
M = 0.04 scales almost as N' 1 (Fig. 4). On the 
other hand, when the loss cone is full, (j>l, the 
decay rate is nearly independent of the relaxation 
time, since (by assumption) scattering occurs so 
quickly that orbits are never depleted. The much 
weaker dependence of s on N for M — 0.01 (Fig. 4) 
is qualitatively consistent with this prediction. 

We tested this model more carefully, by com- 
puting the expected hardening rate of the binary, 
under the assumption that the supply of stars was 
limited by the rate of loss cone repopulation. The 
flux of stars into the binary's loss cone (mass per 
unit time) in a steady-state spherical galaxy is 
given approximately by 



T{E)dE « Aw 2 Jf c (E)q(E) 



f(E) 
In R- 1 



dE 



(11) 



(Paper I), where f(E) is the phase-space mass 
density of stars at energy E, and Rq(E) is the ef- 
fective, dimcnsionless size of the loss cone as seen 
by a star of energy E; an approximate expression 
for Rq(E) is given by Cohn & Kulsrud (1978). If 
all of the stars scattered into the loss cone are as- 
sumed to interact instantaneously with the binary, 
its decay rate is roughly 



s = 



1 

dt V a 



2(g) 
aM 



I 



T(E)dE (12) 



where (C) is the average value of the dimension- 
less energy change during a single star-binary en- 
counter, C = (M/2m*)(AE/E) (Hills 1983). In 
the limit a <C ah, (C) is constant and roughly 
equal to one. While equation (12) with constant 
(C) is unlikely to give a good description of the 
early evolution of the binary, when a ^ ah, it 
should reproduce the late-time hardening rate of 
the binary fairly well. Figure 4 shows that this is 
indeed the case; we set r t = a and (C) — 1.25, 
and evaluated s from equation (12) at the final 



value of a reached in each integration. Given the 
simplicity of the model, the agreement with the 
iV-body hardening rates is remarkably good. De- 
viations between theory and simulation appear to 
be greatest at largest N; this may be due to the 
fact that the binary separation does not fall far 
below ah in these integrations. 

We conclude that the binary hardening rates 
observed in the iV-body simulations - including 
their dependence on M and N - are consistent 
with the predictions of loss cone repopulation the- 
ory. In particular, the integrations with M = 
0.04 exhibit the ~ iV -1 scaling of hardening rate 
with particle number that is characteristic of an 
"empty," diffusively-repopulated loss cone. 

Brownian motion of the binary might also affect 
its decay rate, by enhancing the diffusion of stars 
into the binary's loss cone. In Paper I, a simple 
model was presented for this process. Equation 
(85) from that paper gives the diffusion time (i.e. 
the time for the binary's loss cone to be repopu- 
lated) due to Brownian motion as 



tBr 



(Ka)(2E) 7 / 2 
6ttG 2 {M 1 +M 2 ) 2 A 2 ' 



(13) 



Here Ka is defined as above, as the distance 
from the binary's center within which the grav- 
itational slingshot is effective; E is the orbital 
energy, assuming a Keplerian potential around 
the black hole (the potential from the other stars 
was ignored in this model); and A 2 is the mean 
square acceleration experienced by the binary due 
to gravitational perturbations from stars. Adopt- 
ing equation (76) from Paper I for A, and set- 
ting E = a 2 , the orbital energy of a star near the 
binary's gravitational influence radius, we find in 
A~-body units 



10 2 N 6 M- 2 a-2 



(14) 



where N 6 = N/10 6 , M_ 2 = (M 1 + M 2 )/0.01, and 
a_ 2 = a/0.01. 

This time scale should be compared with the 
time scale for stars to diffuse into the binary's 
loss cone due to star-star encounters, or tdiff ~ 
P(E)/q(E) with P(E) the orbital period. Evalu- 
ating q and P at the binary's gravitational influ- 
ence radius, E w Eh = o 2 , we find 



tdiff « {0.59, 0.15} q(E h y 



(15) 
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for M = {0.04, 0.01}. In the simulations with the 
smaller binary, Figure 8 shows that q is always 
greater than ~ 1, implying tdiff ^ 0.15. This 
is much shorter than tsrown for any (N, a) con- 
sidered here, from which it follows that Brownian 
motion had almost no effect on loss-cone repop- 
ulation. In effect, the smaller binary's loss cone 
was maintained in such a full state by star-star 
encounters that the additional influence of the bi- 
nary's motion was negligible. In the case of the 
more massive binary, we have 

J«£f- n (itfqNea-t)- 1 . (16) 

t Brown 

Figures 2 and 8 suggest that this ratio might have 
approached unity at late times in the simulations 
with smallest N, but generally was much smaller 
than one, again suggesting that the influence of 
Brownian motion on the binary's evolution was 
small. 

5. Comparison with Other Work 

Our initial conditions were chosen to be iden- 
tical with those of Chatterjee, Hcrnquist & Loeb 
(2003) (CHL), who used a rather different TV-body 
algorithm to follow the evolution of massive bi- 
naries at the centers of Plummer-model galaxies. 
These authors present detailed results for only 
two integrations, with M\ = M% = 0.00125 and 
Mi — M 2 = 0.02; the latter values match our 
more massive binary, and Figure 2 in CHL shows 
the evolution of 1/a, e and p(0) for an integration 
with N = 0.2 x 10 6 . The evolution found by CHL 
for e and p(0) appears consistent with what we 
find here for the same values of M and N. How- 
ever our results for the binary hardening rate differ 
from theirs, in two respects. (1) CHL observed bi- 
nary hardening rates s = (d/dt)(l/a) that were 
significantly time-dependent, with 1/a sometimes 
increasing as weakly as ~ t 5 , rather than the lin- 
ear dependence observed here. For instance, their 
Fig. 2 shows a gradually falling decay rate in the 
case (Mi = M 2 = 0.02, N = 0.2 x 10 6 ), with a 
value of 1/a at t = 250 of ~ 360, compared with 
our value of - 500. (2) While CHL also found that 
hardening rates tended to fall with increasing N, 
they state that "s falls systematically until there 
are roughly 200, 000 - 400, 000 stars, when it sta- 
bilizes to a particular value, s ." Unfortunately, 
CHL provided no plots of s(t) or s(N) against 



which this result could be verified; however they 
did quote some particular values for so, e.g. so ~ 8 
when Mi = Mi = 0.005. This is consistent with 
the value of s that we measure for large N in our 
integrations with the same binary mass (Fig. 4), 
but we observe no tendency for the hardening rate 
to "stabilize" at large N, either for this binary or 
for the more massive one that we integrated. 

CHL compared the results of their integrations 
with the predictions of local theory, via an equa- 
tion similar to our equation (9). Since that equa- 
tion predicts no dependence of hardening rates on 
M or N, CHL invoked the mass dependence of the 
binary's Brownian motion to explain their results. 
CHL suggested (without presenting a quantitative 
theory) that the larger wandering radius in inte- 
grations with smaller M and N would translate 
into larger rates of interaction between the binary 
and the stars. The "stabilization" of the hard- 
ening rate which they observed at large N (and 
which we do not observe) was attributed to a kind 
of Brownian- motion- mediated feedback process, in 
which the binary maintains a constant supply rate 
by modulating the local density of stars. However 
no supporting evidence for this model was pre- 
sented; for instance, it was not demonstrated that 
the central density was actually regulated by the 
binary in their TV-body integrations, or that the 
amplitude of the Brownian wandering increased 
with N in the manner postulated. 

Decay of binary black holes in Plummer-model 
galaxies was also investigated by Hemsendorf et 
al. (2002) (HSS). Those authors placed two equal- 
mass particles, Mi = M2 = 0.01, at the center of 
the galaxy and carried out integrations with par- 
ticle numbers N = (0.033, 0.065, 0.131) x 10 6 . The 
black hole particles were placed on initially eccen- 
tric orbits at radii ±0.5, farther from the center 
than in our or in CHL's simulations. Furthermore 
the models were integrated only until t = 60, and 
in the case of the largest N, only until t « 40. 
Comparison with our Figure 3 suggests that the 
hardening rate of the binary at t « 60 might 
be rather different from its late-time value. HSS 
found a mean hardening rate for all of their inte- 
grations of s w 6 — 7 (their Fig. 3), quite con- 
sistent with the early hardening rates found here 
(Fig. 3), especially if one takes into account that 
their binary mass lies midway between our two. 
HSS observed only a slight dependence of decay 
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rate on TV, but this too is consistent with what we 
find at early times (Fig. 3). HSS observed a much 
stronger evolution of the binary's eccentricity, but 
this may be due to their large initial eccentricity, 
e w 0.8 — 0.9. Eccentric binaries are expected to 
evolve in the direction of increasing eccentricity 
(Mikkola & Valtoncn 1992; Quinlan 1996). 

A number of other iV-body studies of binary 
evolution have been carried out using galaxy mod- 
els with large cores (Ebisuzaki ct al. 1991; Makino 
ct al. 1993; Governato et al. 1994), though most 
of these used too few particles to show a depen- 
dence of hardening rate on N. One exception was 
the recent study by Funato & Makino (2003), who 
considered binaries in King (1966)-model galaxies 
with particle numbers as high as 1.0 x 10 6 . These 
authors also failed to reproduce the stabilization 
of the binary hardening rate at large ./V found by 
CHL. 

6. Implications for Binary Evolution in 
Real Galaxies 

The strong TV-dependence of the binary hard- 
ening rate in our simulations makes it dangerous 
to extrapolate our results to real galaxies. Such 
an extrapolation may nevertheless be justified in 
the case of our more massive binary, which ap- 
peared (Fig. 8) to be in the "empty loss cone" 
regime. Massive binaries in real galaxies are also 
expected to be in this regime (Paper I) , at least in 
the absence of any additional mechanisms for loss- 
cone repopulation, and the ./V-dependence of the 
hardening rate observed here for the more massive 
binary might therefore be valid for much larger N. 

We begin by noting that the degree of evolu- 
tion, as measured by the ratio ah/a(tf) between 
initial and final binary separations, was equal to 
- (21,14,8.8,6.2) for N = (0.05,0.1,0.2,0.4) x 
10 6 in the simulations with the more massive bi- 
nary. If the binary's actual mass was 1O 8 M , the 
M — a relation (Ferrarese & Ford 2005) implies 
a galaxy velocity dispersion of a « 180 km s _1 , 
so ah w 1.7 pc and the final scaled separations be- 
come (0.025, 0.12, 0.19, 0.27) pc. Coalescence of an 
equal-mass binary in a Hubble time due to grav- 
itational wave emission requires a^/a J> 20 (Mer- 
ritt & Milosavljevic 2004), so the integration with 
smallest N can be said to have just reached the 
gravitational radiation regime. 



But the implied degree of hardening in a real 
galaxy is much less than this if we take into ac- 
count the likely scaling of the hardening rate with 
N. Ignoring possible changes in the central den- 
sity of the host galaxy (such changes were small 
in the integrations reported here), equation (9a) 
implies 



r d 








dt 





-1 



ss 2.0 x 10 4 [Gp(0)p 1/2 (^) N\f%) 

or 

with T cr = (Gp(0)) -1 / 2 the galaxy crossing time. 
Equation (18) suggests that binary hardening 
times in real galaxies would be very long, of order 
10 7 crossing times and much greater than the age 
of the universe. In effect, the binary would "stall." 

A very different conclusion was reached by 
Chatterjee et al. (2003), who assumed that the bi- 
nary hardening rates observed by them at N 
0.2 — 0.4 x 10 6 would remain constant even for 
much larger values of N. 

Of course our prediction is based on a rather 
unphysical (too homogeneous) galaxy model, an 
unrealistic (too large) choice for the mass of the 
binary, and unlikely initial conditions (with the 
two black holes deposited symmetrically about the 
center of the galaxy, rather than arriving there as 
a result of a galactic merger). Binary hardening 
rates in galaxies with steeply-rising central den- 
sities might be significantly larger than implied 
by equation (18). Upcoming papers in this se- 
ries will attempt to correct these deficiencies with 
more realistic galaxy models and more physically- 
motivated initial conditions. 

7. Conclusions 

The evolution of a massive binary at the cen- 
ter of a galaxy with a large, constant-density core 
was followed using a high-accuracy iV-body code 
and a special-purpose computer cluster. Two val- 
ues for the binary mass were considered: (Mi + 
M 2 )/M ga i = (0.04,0.01). The binary's hardening 
rate exhibited a clear TV-dependence, consistent 
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with the predictions of collisional loss-cone repop- 
ulation theory. In the simulations with M1 + M2 = 
0.04M g ai and large particle number, TV ^ 0.2 x 10 6 , 
the binary appeared to be nearly in the "empty 
loss cone" regime believed to characterize real 
galaxies (in the absence of other mechanisms for 
loss-cone refilling): the hardening rate decreased 
almost as steeply as TV" 1 . We observed no indica- 
tion that the binary hardening rate "stabilizes" at 
particle numbers greater than TV « 0.2 x 10 6 , or 
that the binary's evolution is strongly affected by 
its Brownian motion, as claimed in a recent study 
based on a more approximate TV-body treatment 
(Chatterjee, Hernquist & Loeb 2003). 

Our results demonstrate the feasibility of TV- 
body simulations of binary black hole evolution 
that mimic the behavior expected in systems of 
much larger TV. Even larger particle numbers, 
N <; 10 6 , will be required to correctly reproduce 
the evolution of binary black holes in more real- 
istic galaxy models with higher central densities 
and less massive binaries. Direct-summation TV- 
body simulations with such high particle numbers 
are now becoming feasible with the advent of par- 
allel, special-purpose computers like gravitySimu- 
lator. Future papers in this series will present the 
results of such simulations. 
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